{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** This example has been significantly expanded and enhanced. The new, recommended version is located [here](solutions/adv3_panorama-stitching-solution.ipynb). We retain this version intact as it was the exact example used in [the scikit-image paper](https://peerj.com/articles/453/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Panorama stitching"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage import io, transform\n",
    "from skimage.color import rgb2gray\n",
    "from skdemo import imshow_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ic = io.ImageCollection('../images/pano/DFM_*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``ImageCollection`` class provides an easy way of\n",
    "loading and representing multiple images.  Images are not\n",
    "read from disk until accessed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow_all(ic[0], ic[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Credit: Photographs taken in Petra, Jordan by François Malan<br/>\n",
    "License: CC-BY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "image0 = rgb2gray(ic[0][:, 500:500+1987, :])\n",
    "image1 = rgb2gray(ic[1][:, 500:500+1987, :])\n",
    "\n",
    "image0 = transform.rescale(image0, 0.25)\n",
    "image1 = transform.rescale(image1, 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow_all(image0, image1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this demo, we estimate a projective transformation\n",
    "that relates the two images.  Since the outer\n",
    "parts of these photographs do not comform well to such\n",
    "a model, we select only the central parts.  To\n",
    "further speed up the demonstration, images are downscaled\n",
    "to 25% of their original size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Feature detection and matching\n",
    "\n",
    "\"Oriented FAST and rotated BRIEF\" features are detected in both images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.feature import ORB, match_descriptors\n",
    "\n",
    "orb = ORB(n_keypoints=1000, fast_threshold=0.05)\n",
    "\n",
    "orb.detect_and_extract(image0)\n",
    "keypoints1 = orb.keypoints\n",
    "descriptors1 = orb.descriptors\n",
    "\n",
    "orb.detect_and_extract(image1)\n",
    "keypoints2 = orb.keypoints\n",
    "descriptors2 = orb.descriptors\n",
    "\n",
    "matches12 = match_descriptors(descriptors1, descriptors2, cross_check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.feature import plot_matches\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 10))\n",
    "plot_matches(ax, image0, image1, keypoints1, keypoints2, matches12)\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each feature yields a binary descriptor; those are used to find\n",
    "the putative matches shown.  Many false matches are observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Transform estimation\n",
    "\n",
    "To filter matches, we apply RANdom SAMple Consensus (RANSAC),\n",
    "a common method of rejecting outliers.  This iterative process\n",
    "estimates transformation models based on\n",
    "randomly chosen subsets of matches, finally selecting the\n",
    "model which corresponds best with the majority of matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.transform import ProjectiveTransform\n",
    "from skimage.measure import ransac\n",
    "from skimage.feature import plot_matches\n",
    "\n",
    "# Select keypoints from the source (image to be registered)\n",
    "# and target (reference image)\n",
    "src = keypoints2[matches12[:, 1]][:, ::-1]\n",
    "dst = keypoints1[matches12[:, 0]][:, ::-1]\n",
    "\n",
    "model_robust, inliers = ransac((src, dst), ProjectiveTransform,\n",
    "                               min_samples=4, residual_threshold=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 15))\n",
    "plot_matches(ax, image0, image1, keypoints1, keypoints2, matches12[inliers])\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how most of the false matches have now been rejected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Warping\n",
    "\n",
    "Next, we want to produce the panorama itself.  The first\n",
    "step is to find the shape of the output image, by taking\n",
    "considering the extents of all warped images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.transform import SimilarityTransform\n",
    "\n",
    "r, c = image1.shape[:2]\n",
    "\n",
    "# Note that transformations take coordinates in (x, y) format,\n",
    "# not (row, column), in order to be consistent with most literature\n",
    "corners = np.array([[0, 0],\n",
    "                    [0, r],\n",
    "                    [c, 0],\n",
    "                    [c, r]])\n",
    "\n",
    "# Warp the image corners to their new positions\n",
    "warped_corners = model_robust(corners)\n",
    "\n",
    "# Find the extents of both the reference image and the warped\n",
    "# target image\n",
    "all_corners = np.vstack((warped_corners, corners))\n",
    "\n",
    "corner_min = np.min(all_corners, axis=0)\n",
    "corner_max = np.max(all_corners, axis=0)\n",
    "\n",
    "output_shape = (corner_max - corner_min)\n",
    "output_shape = np.ceil(output_shape[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warp the images according to the estimated transformation model.\n",
    "Values outside the input images are set to -1 to distinguish the\n",
    "\"background\".\n",
    "\n",
    "A shift is added to make sure that both images are visible in their\n",
    "entirety.  Note that ``warp`` takes the inverse mapping\n",
    "as an input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.color import gray2rgb\n",
    "from skimage.exposure import rescale_intensity\n",
    "from skimage.transform import warp\n",
    "\n",
    "offset = SimilarityTransform(translation=-corner_min)\n",
    "\n",
    "image0_ = warp(image0, offset.inverse,\n",
    "               output_shape=output_shape, cval=-1)\n",
    "\n",
    "image1_ = warp(image1, (model_robust + offset).inverse,\n",
    "               output_shape=output_shape, cval=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alpha channel is now added to the warped images\n",
    "before they are merged together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def add_alpha(image, background=-1):\n",
    "    \"\"\"Add an alpha layer to the image.\n",
    "    \n",
    "    The alpha layer is set to 1 for foreground and 0 for background.\n",
    "    \"\"\"\n",
    "    return np.dstack((gray2rgb(image), (image != background)))\n",
    "\n",
    "image0_alpha = add_alpha(image0_)\n",
    "image1_alpha = add_alpha(image1_)\n",
    "\n",
    "merged = (image0_alpha + image1_alpha)\n",
    "alpha = merged[..., 3]\n",
    "\n",
    "# The summed alpha layers give us an indication of how many\n",
    "# images were combined to make up each pixel.  Divide by the\n",
    "# number of images to get an average.\n",
    "merged /= np.maximum(alpha, 1)[..., np.newaxis]\n",
    "merged = merged[..., :3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow_all(image0_alpha, image1_alpha, merged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, while the columns are well aligned, the color\n",
    "intensity is not well matched between images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Blending\n",
    "\n",
    "To blend images smoothly we make use of the open source package\n",
    "[Enblend](http://enblend.sf.net), which in turn employs multi-resolution splines and \n",
    "Laplacian pyramids [1, 2].\n",
    "\n",
    "[1]\tP. Burt and E. Adelson. [\"A Multiresolution Spline With Application to Image Mosaics\"](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.56.690). ACM Transactions on Graphics, Vol. 2, No. 4, October 1983. Pg. 217-236.\n",
    "[2]\tP. Burt and E. Adelson. [\"The Laplacian Pyramid as a Compact Image Code\"](https://doi.org/10.1109/TCOM.1983.1095851). IEEE Transactions on Communications, April 1983."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imsave('/tmp/frame0.tif', image0_alpha)\n",
    "plt.imsave('/tmp/frame1.tif', image1_alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "enblend /tmp/frame*.tif -o /tmp/pano.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pano = io.imread('/tmp/pano.tif')\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(pano)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "<div style=\"height: 400px;\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%reload_ext load_style\n",
    "%load_style ../themes/tutorial.css"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
